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ABSTRACT 


Atmospheric  turbulence  will  induce  phase  and  amplitude  fluctuations  in 
propagating  electromagnetic  waves,  such  as  a  laser  beam.  Adaptive  optical  systems 
attempt  to  compensate  for  these  distortions.  The  Strehl  ratio  is  a  measure  of  the  peak,  on- 
axis  intensity  after  propagation  through  turbulence  divided  by  the  peak  irradiance  for 
vacuum  propagation. 

This  thesis  investigated  the  probability  distribution  of  the  Strehl  ratio  of  a  perfect, 
phase-only,  adaptive  optical  system  as  a  function  of  the  atmospheric  coherence  length,  r0, 
divided  by  the  actuator  spacing,  d.  Using  an  efficient  Fourier  algorithm  and  28  work- 
stations running  in  parallel,  over  850  million  computer  simulations  were  performed  for  25 
different  d/r0  ratios  in  order  to  produce  a  histogram  of  the  Strehl  ratio  probability 
distribution.  The  results  show  that  the  Strehl  ratio  follows  a  log-normal  probability 
distribution  even  for  very  small  probabilities. 

A  second  set  of  computer  simulations  introduced  intensity  scintillation  by 
including  the  log-amplitude  variance  parameter,  g(  .  Much  faster,  state-of-the-art 
computer  workstations  enabled  over  two  billion  realizations  on  1 8  machines  running  in 
parallel  for  comparable  time  periods.  The  trends  of  these  results  are  more  complex  and 
will  require  further  research  and  deeper  investigation. 
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I.         INTRODUCTION 

A  laser  beam  propagating  through  the  atmosphere  will  experience  degradation  due 
to  atmospheric  turbulence.  Temperature  differences  associated  with  the  turbulence 
generate  changes  in  the  air  density,  and  subsequently  its  index  of  refraction,  with  the 
result  that  parts  of  the  laser  beam's  spherical  wavefront  are  slowed  by  differing  amounts. 
The  field  of  adaptive  optics  uses  a  deformable  mirror  or  a  similar  device  to  compensate, 
or  correct,  for  the  distortion  of  the  light,  effectively  reversing  the  turbulent  effects  and 
restoring  the  spherical  shape  of  the  wavefront. 

There  are  three  important  parameters  of  significance  to  this  problem.  The  first  is 
the  coherence  length,  or  r0,  which  is  a  measurable  parameter  that  indicates  the  severity  of 
atmospheric  turbulence  induced  wavefront  distortion,  and  it  is  extremely  variable 
depending  on  atmospheric  conditions  at  a  particular  site.  The  second  is  the  log-amplitude 
variance,  ce  ,  which  is  a  measure  of  the  intensity  variation  of  the  light  perturbed  by 
turbulence. 

The  third  parameter  is  the  Strehl  ratio,  S,  which  is  a  measure  of  a  light  beam's 
intensity  at  the  nominal  focus,  relative  to  what  it  would  be  if  the  focusing  process  were 
diffraction  limited.  The  best  one  can  hope  for  is  to  achieve  diffraction-limited  intensity, 
which  can  be  approached  with  very  good  adaptive  optics  correction  of  atmospheric 
turbulence  effects  -  but  which  cannot  be  reached. 


Pertinent  to  this  thesis  is  the  relationship  between  r0  and  <je  and  the  randomly 
varying  values  of  the  Strehl  ratio  of  a  laser  beam  degraded  by  atmospheric  turbulence  and 
(partially)  compensated  by  adaptive  optics.  The  nominal  (or  average)  value  of  the  Strehl 
ratio  will  be  degraded  as  r0  gets  smaller  and  as  <je  gets  larger.  However  the  value  of  the 
Strehl  ratio  at  any  instant  is  a  random  variable  and  there  is  no  exact  relationship  between 
that  value  and  the  values  of  r0  and  <je  .  In  fact,  for  given  values  of  r0  and  <je  there  will  be 
many  different  patterns  of  turbulence  effects  on  the  light  beam,  each  of  which  will  yield 
its  own  Strehl  ratio.  A  study  of  the  nature  of  the  distribution  of  the  random  values  of  the 
Strehl  ratio  and  the  relation  of  the  distribution  to  the  values  r0  and  oe  is  the  subject  of 
interest  here. 

The  problem  of  assessing  the  magnitude  and  significance  of  the  Strehl  ratio 
degradation  experimentally  is  made  very  difficult  by  the  time  and  effort  required 
collecting  these  measurements  in  actual  site  situations.  Not  only  is  the  required 
equipment  complex  and  expensive,  but  also  conditions  change  so  rapidly  with  time  and 
one  can  never  "select"  particular  weather  conditions  to  accumulate  data  measurements. 
As  a  consequence  a  very  great  deal  of  time  in  the  field  would  be  required  to  implement 
such  an  approach. 

Fortunately,  due  to  exponential  advances  in  computer  technology,  computer 
simulation  modeling  is  becoming  increasingly  viable  for  all  fields  of  endeavor.  This 
thesis  will  simulate  the  propagation  of  a  laser  beam  through  the  atmosphere,  with 
selectable    values    of  r0    and    oe     and    then    evaluate    the    subsequent    Strehl    ratio 


measurements.  By  combining  the  computer  resources  across  the  Naval  Postgraduate 
School,  I  have  been  able  to  run  as  many  as  twenty-eight  Monte  Carlo  simulations  in 
parallel  and  generate  very  large  sets  of  the  Strehl  ratio  values  from  which  the  probability 
distributions  can  be  calculated  with  statistically  significant  results  at  very  low  probability 
levels. 


II.       BACKGROUND 

A.        OVERVIEW 

The  atmosphere  has  multiple  affects  on  the  propagation  of  a  laser  beam. 
Molecular  and  aerosol  scattering  and  absorption  will  attenuate  the  beam.  High  intensity 
laser  beams  can  heat  the  atmosphere  as  it  travels  -  with  consequent  effects  on  the 
propagation  of  the  laser  beam,  which  effects  are  known  as  thermal  blooming.  In  addition, 
atmospheric  turbulence  will  distort  the  phase  of  a  beam  and  make  its  intensity  non- 
uniform as  a  result  of  random  variation  of  the  atmosphere's  index  of  refraction  and  multi- 
path  interference  induced  by  the  turbulence.  As  this  thesis  will  only  address  lower  power 
laser  beams  not  significantly  affected  by  thermal  blooming,  atmospheric  turbulence  will 
be  the  only  factor  considered  (Greenwood,  1992,  pp. 3-7). 

Figures  (1)  and  (2)  on  the  following  pages  show  the  effect  of  the  atmosphere  on 
laser  propagation.  These  figures  plot  the  intensity  distribution  as  a  hot-cold  differential. 
In  Figure  (1),  the  intensity  appears  as  a  hot-cold  color-scale  for  a  Gaussian  beam  after 
long  distance  propagation  through  a  vacuum.  The  diameter  of  the  Gaussian  waist  will  be 

proportional  in  this  case  to  the  wavelength  divided  by  the  aperture  diameter   (w=yL). 

However,  after  passing  through  a  turbulent  medium  as  Figure  (2)  illustrates,  the  resulting 
beam  pattern  includes  spreading,  intensity  and  phase  variations. 


Figure  1.  Laser  Intensity  Distribution  through  a  Vacuum 


B.         COHERENCE  LENGTH 

At  the  heart  of  describing  atmospheric  turbulence  is  the  index  of  refraction 
structure  parameter.    Denoted  byC„2,  it  relates  the  mean  square  difference  in  index  of 

refraction  at  two  arbitrary  points  to  the  two-thirds  power  of  the  distance  between  the 
points.  It  is  defined  by  the  following  equation:  (Tatarski,  1961) 


C    = 


[n2-nxy 


.2/3 


(1) 


Here,  r  denotes  the  distance  between  the  two  points  and  the  angled  brackets  indicate  an 


Figure  2.  Laser  Intensity  Fluctuations  through  a  Turbulent  Medium 
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ensemble  average  value.    The  r     dependence  arises  from  Kolmogorov  statistics  of  the 
turbulence. 


In  practice,  it  is  usually  more  convenient  to  measure  C„2  indirectly,  by  measuring 
the  analogous  temperature  structure  parameter,  C2T  and  relating  back  to  C„2  according  to 
the  relation 


Cl  = 


r  p\2 


79-10" 


V 
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c 


r> 


(2) 


where  P  is  the  atmospheric  pressure  in  Pascal  and  T  is  the  atmospheric  temperature  in 


Kelvin.   Balloon-launched  devices  with  a  pair  of  fast-response  temperature  probes 
attached  are  commonly  used  to  measure  C\ ,  and  from  this,  C„2  values  are  determined. 

Once  C„2  is  known  along  the  path  of  beam  propagation,  the  coherence  length,  or  r0,  can 

be  calculated  through  the  following  equation  for  a  ground  to  space  propagation:    (Fried, 
1966) 


r„  =2.1 


z 

1.46&2sec6>  \C2n{z)dz 


(3) 


In  this  equation,  k  =  2nl  X  is  the  wavenumber  and  #is  the  zenith  angle  of  the  path  whose 
length  is  z.  The  coherence  length  is  a  measure  of  the  electric  field  auto-correlation  length 
and  indicates  the  severity  of  atmospheric  turbulence.  It  is  extremely  variable  depending 
on  atmospheric  conditions  at  a  particular  site. 

The  significance  of  the  coherence  length  arises  with  an  example  from  astronomy 
(for  light  propagating  down  through  the  atmosphere),  although  it  applies  exactly  the  same 
for  laser  beam  propagation  up  through  the  atmosphere,  through  the  principle  of 
reciprocity.  When  considering  the  resolution  of  a  telescope,  it  is  a  well-understood  fact 
that  in  the  absence  of  turbulence,  the  angular  separation  ((psep)  of  two  stars  that  can  just  be 
resolved  is  proportional  to  the  wavelength  divided  by  the  aperture  diameter  (Born  & 
Wolf,  1970,  p.415).  For  conventional  telescopes,  r0  is  the  diameter  of  the  largest  aperture 
that  can  be  used  without  turbulence  significantly  degrading  the  image.   The  effect  is  that 

X        X 

<p     ex >  — .    As  the  turbulence  gets  stronger,  the  values  for  r0  are  driven  smaller  and 

D       ro 


the  resolution  limit/separation  (psep  will  get  larger/coarser.  For  earthbound  observatories 
operating  at  visible  wavelengths,  r0  typically  ranges  between  five  and  fifteen  centimeters 
(Hardy,  1994,  p.61). 

In  addition  to  the  spreading  of  an  image  spot  by  random  phase  variations, 
turbulence  also  induces  variations  in  the  intensity  pattern.  This  is  a  process  perhaps  most 
simply  characterized  as  "random  apodization."  Both  phase  and  intensity  variations  effect 
the  size  of  the  image  spot  and  the  system's  Strehl  ratio. 

C.         STREHL  RATIO 

There  have  been  literally  hundreds  of  volumes  written  on  the  propagation  of  light 
and  its  mathematical  characteristics,  so  only  a  brief  summary  is  provided  here,  pertinent 
to  the  problems  addressed  in  this  thesis.  Figure  (3)  below  represents  a  coordinate  system 
for  discussion  (Tyson,  1991,  p. 6).  As  a  means  to  simplify  terms  in  the  exponent  of  forth- 


Figure  3.  Reference  Coordinates  for  Electromagnetic  Propagation 


coming  expressions,  the  following  normalized  coordinates  in  the  focal  plane  are  used: 


u  = 
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z, 


(4)        and 
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For  a  beam  of  coherent  light  at  wavelength  A,  and  variable  amplitude  function,  Ae  Cip,&) , 
where  £  =  rj  +  i(j)  can  be  complex,  the  intensity  of  light  at  a  point  P  on  the  image  plane  a 
distance  z  away  is  given  by:  (Born  &  Wolf,  1970,  p.461) 


/(/>)  = 


^a2 
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\\ 

0    0 


2     ]pdpd6 
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The  notation  used  here  is  illustrated  in  Figure  (3),  where  the  variables  (p,0,z)  constitute 
a  cylindrical  coordinate  system.  R  is  the  slant  range  from  the  center  of  the  pupil  to  point 
P,  <p  represents  the  deviation  in-phase  from  a  perfect  sphere  about  the  origin  of  the  focal 
plane,  and  rj  represents  the  intensity  deviations. 

From  the  standpoint  of  adaptive  optics,  and  pertinent  to  this  thesis,  the  most 
important  quantity  in  the  above  expression  is  the  cj)  term,  commonly  referred  to  as  the 
wavefront  phase  perturbation.  If  no  aberrations  are  present  (diffraction-limited  case),  then 
the  intensity  is  a  maximum  on-axis  (r=0),  also  called  the  Gaussian  image  point:  (Born  & 
Wolf,  1970,  pp.46 1-2) 


Vo=^=o=^ 


(  Aa1^ 


KXR 


(7) 


The  Strehl  ratio,  S,  is  a  measure  of  a  beam's  intensity  at  focus  (or  in  the  far- field) 
relative  to  what  it  would  be  if  there  were  no  optical  system  induced  or  turbulence  induced 
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aberrations.  It  is  defined  as  the  ratio  between  the  intensity  at  focus  or  on-axis  in  the  far- 
field  of  an  aberrated  beam  and  the  intensity  of  an  unaberrated  beam  at  the  same  position. 
The  remaining  aberrations  are  represented  by  £  and  the  Strehl  ratio  can  be  expressed  as: 
(Born  &  Wolf,  1970,  pp.462) 


_  C\\  _      ^    Origin  )  _ 


S(v  =  0)  = 


1  In 


\\eQ  pdpdO 


0   0 


(8) 


For  phase  only  adaptive  optics,  intensity  variations  are  ignored  (because  they 
cannot  be  corrected),  and  rj  is  often  set  to  zero.  If  both  the  intensity  and  phase  effects  of 
the  beam  at  the  pupil  are  unaberrated  (or  (fa  =  zero),  the  Strehl  ratio  will  reduce  to  unity, 
S  =  1 ,  which  is  to  say  that  the  intensity  at  the  focal  plane  is  diffraction-limited.  This  is 
the  unattainable  goal  of  perfect  adaptive  optics;  to  correct  the  wavefront  in  such  a  manner 
that  (f>p  =  0  and  the  Strehl  ratio  is  unity. 

D.         COMPUTER  SIMULATION  OF  BEAM  PROPAGATION 

A  word  needs  to  be  said  concerning  the  implementation  of  beam  propagation 
using  computer  simulation.  There  are  essentially  two  methods.  The  first  method  is  most 
intuitive  to  the  casual  observer,  by  defining  a  beam  profile  and  actually  propagating  it  by 
using  Fourier  transform  matrices  as  the  Optical  Transfer  Functions  (OTF's).  One  can 
create  the  random  phase  screen  (representing  the  turbulence),  then  take  the  Fourier 
transform  and  multiply  by  the  input  beam  pattern.  The  propagation  distance  is  important, 
because  the  underlying  principle  is  Huygens-Fresnel  propagation.     This  method  of 
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simulation  can  be  extremely  accurate,  although  multiple  passes  are  frequently  necessary 
to  represent  long  propagation  lengths  through  different  turbulent  areas  (Davis,  1 994,  pp. 
34-36).  One  other  important  point  is  the  fact  that  both  amplitude  and  phase  effects  from 
the  turbulence  occur  with  each  propagation  step.  For  the  second  method,  the  amplitude 
effects  must  be  characterized  separately  from  the  phase  effects,  and  then  combined  at  a 
later  time  by  use  of  convolution  (Fried,  Sep  1997,  p.ii). 

The  second  method  is  based  on  reciprocity  in  optical  propagation.  It  allows  us, 
instead  of  considering  propagation  of  a  laser  beam  to  the  far-field  through  multiple  layers 
of  turbulence  -  instead  to  consider  a  distorted  laser  beam  arriving  at  an  aperture  after 
propagation  through  turbulence  and  then  to  simply  consider  the  propagation  of  that  laser 
beam,  in  one  step,  to  focus.  Use  of  this  second  method  allows  us  to  make  use  of  our 
knowledge  of  the  statistics  of  the  beam's  distortion  after  propagation  through  the 
distributed  turbulence  to  avoid  the  multi-step  calculation  of  propagation  through 
distributed  turbulence.  Ignoring  intensity  variations,  which  we  do  in  this  section  (though 
not  in  the  next),  we  simulate  the  phase  variations  of  the  laser  beam  after  it  has  propagated 
through  the  distributed  turbulence  by  a  single  phase  screen.  (To  have  simulated 
propagation  through  the  distributed  turbulence,  we  would  have  needed  a  number  of  such 
-  weaker  -  phase  screens.) 

The  simulation  of  the  phase  screen  requires  knowledge  of  the  power  spectral 
density  of  the  medium  through  which  the  beam  will  propagate.  This  context  requires  an 
accurate  power  spectral  density  for  the  phase  screens.  Fortunately,  much  research  has 
been  completed  in  this  area  by  pioneers  in  the  field,  such  as  Tatarski,  Kolmogorov,  and 
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Fried.  Without  showing  the  full  development,  it  can  be  shown  that  since  the  phase 
structure  function  follows  a  five-thirds  power  law,  then  the  power  spectral  density,  (£>p , 
has  the  form:  (Fried,  Sep  1997,  p.2) 

0>  (K)  =  Ar0-5"\K[n'\  (9) 


where 
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A=-7T-*13 

3 


4F5 


5/6       ,       v       ,      s-i 


r 


r 


*  0.045791 17422.  (10) 


voy 


This  equation  for  the  power  spectral  density  provides  a  means  for  generating  phase 
screens  that  contain  the  proper  atmospheric  phase  statistics.  This  formulation  allows  easy 
scalability  to  different  turbulence  conditions  through  the  r0-term. 

Section  III  of  this  thesis  includes  the  effects  of  the  turbulence  along  the 
propagation  path  in  terms  of  the  wavefront  distortion  at  the  optics  aperture,  without  the 
effects  of  intensity  scintillation  included.  The  value  of  r0  characterizes  the  wavefront 
distortion  and  the  phase  structure  has  a  five-thirds  power-law  dependence  (Fried, 
Sep.  1997,  p.l).  The  efficiency  of  this  method  is  greater  than  the  first  method  described 
above,  because  the  entire  propagation  can  be  represented  by  one  phase  screen,  requiring 
only  one  Fourier  transform  to  complete  the  simulated  propagation.  Additionally  by 
convolution,  a  two-dimensional  phase  screen  (256  x  256)  contains  separate  realizations. 

The  effects  of  intensity  scintillation  can  be  simulated  separately  and  in  the  same 
manner  as  the  phase  effects,  which  is  to  say  finding  an  expression  to  represent  the  power 
spectral  density  for  the  intensity  scintillation  pattern.    Intensity  scintillation  is  an  effect 
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that  cannot  be  corrected  by  any  conventional  form  of  adaptive  optics  (Fried,  May  1997,  p. 
ii).  This  would  be  extremely  troublesome  except  for  the  fact  that  intensity  scintillation 
effects  in  the  atmosphere  are  much  smaller  than  the  phase  effects  (Greenwood,  1979,  p. 
550).  Hence,  adaptive  optic  systems  are  proving  very  effective  by  just  compensating  for 
the  phase  distortions.  Nevertheless,  the  effects  of  intensity  scintillation  are  non-negligible 
and  will  be  addressed,  later  in  Chapter  IV. 
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III.     STREHL  DEGRADATION  DUE  TO  PHASE  EFFECTS 

A.        OVERVIEW 

These  investigations  utilized  computer  simulation  source  code  written  by  Dr. 
David  Fried  in  the  MATLAB  programming  language,  and  as  such,  required  modifications 
before  applying  to  this  particular  problem.  In  the  course  of  accomplishing  the 
modifications,  a  considerable  amount  of  effort  was  required  to  understand  the  intricate 
programming  techniques  used  to  arrive  at  an  efficient  algorithm.  After  all.  the  ultimate 
task  is  to  produce  computer  simulation  code  that  can  produce  very  large  samples  of 
random  strehl  ratio  measurements  efficiently  in  order  to  investigate  very  small  probability 
events. 

An  example  here  is  appropriate:  Consider  a  turbulence  strength  of  interest  being 
an  average  value  (r0  =  10  cm),  and  one  measurement  of  the  strehl  ratio  from  an  adaptive 
optics  compensated  system  results  to  S  =  0.35.  Is  this  average?  Is  it  extremely  good  and 
luck  was  on  our  side?  Or  is  it  much  worse  than  had  been  hoped,  and  luck  definitely  was 
not  on  our  side?  One  very  good  way  to  answer  this  question  is  to  collect  a  large  number 
of  sample  measurements  of  the  strehl  ratio  when  the  r0  value  is  always  10  cm,  and  then 
look  at  the  average  and  standard  deviation  of  the  sample.  The  larger  the  number  of 
samples,  the  more  confident  one  can  be  in  saying  that,  for  example,  the  0.35  strehl  ratio 
was  indeed  an  average  and  there  is  a  significant  probability  that  the  next  measurement 
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will  be  as  poor  as  0.25,   ...   or  not.  It  is  the  efficient  development  of  very  large  numbers 
of  sample  values  that  is  the  core  of  the  work  reported  here. 


B.         PHASE  SCREEN  GENERATION 

To  generate  a  random  phase  screen,  one  must  first  generate  a  suitable  'random' 
starting  set  of  numbers,  and  then  scale  to  the  desired  statistics.  Remembering  that  the 
spatial  frequency  components  are  continuous,  a  substitution  will  be  made  to  account  for 
discrete  representation  in  the  computer  simulation  environment.  Discrete  spatial 
frequencies  will  be  noted  by  '  in  '  in  place  of  the  continuous  '  k  '  spatial  frequency,  both 
in  units  of  cycles/unit  length.  As  such,  '  wim)'  will  be  defined  to  be  the  N  x  N  array  of 
random  values  -  each  value  with  its  corresponding  one  of  the  N  x  TV  spatial  frequencies. 
In  the  phase  screen  generation  algorithm  that  is  used  here,  the  values  of  w(in)  are  chosen 
from  a  Gaussian  distribution  with  zero  mean  and  unity  variance  -  in  other  words, 
standard  white  noise.  The  values  of  w(in)  will  be  complex:  (Fried,  1996,  p.3) 

w(in)  =  wR  (in)  +  iw,  (in). ,  (11) 

with  wR  (in)  and  w7  (in)  each  being  real  Gaussian  random  variables  with  unity  variance. 

As  stated  previously,  from  equations  (9)  and  (10),  the  power  spectral  density 
produces  a  suitable  random  phase  screen.  That  is  to  say  that  the  phase  screens  generated 
will  have  a  five-thirds  power  law  structure  function.  For  this  to  be  so,  from  equation  (9), 
the  spatial  frequency  components  of  the  power  spectral  density  must  be  proportional  to 
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the  minus  eleven-thirds  power  of  the  magnitude  of  the  spatial  frequency  (Fried,  1 996. 
pp. 2-3).  Accordingly,  we  write 


<f>p(m)cc^\m\~        ,  (12) 

where,  because  the  power  spectral  density  will  be  proportional  to  the  square  of  the 
magnitude,  to  arrive  with  a  minus  eleven-thirds  power  law,  the  square  root  has  been 
applied:  (Fried,  1996,  p.4) 

,  r_v      {a\m\-U'6w{m),     if\m\*0 

<f>p{m)=\    II  v  .     _  .  (13) 

[  0  ,  if\m\  =  0 

Here,  the  a-coefficient  encompasses  the  ability  to  scale  the  'strength'  of  the  phase  screen 

so  that  it  will  correspond  to  the  desired  value  r0,  and  will  be  discussed  in  the  next 

paragraph.   Also,  realizing  that  the  zero  spatial  frequency  component  will  only  affect  the 

average  (or  dc)  value  of  the  simulated  phase  screen  it  has  been  set  to  zero,  avoiding  an 

infinite  result. 

The  value  of  the  or-coefficient  is  developed  and  given  by  Fried  to  be  dependent  on 

the  size,  N,  of  the  array:  (Fried,  1996,  pp.  13- 15) 

N-5n-juN-2 
a  =  fi ,  (14) 

where  //  =  0.98965794  ,  and  J3=  6.8838772,  and  the  values  for  S2(N)  and  S,(N)  are  given 
in  Table  (1). 

The  random  array  <pp (m)  conforming  to  the  power  spectral  density  characteristics 
of  a  random  turbulent  phase  perturbation,  can  be  converted  back  into  the  spatial  domain 
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TV 

Sj(N) 

S2(N) 

16 

2.27148846e-07 

2.26348792e-05 

32 

3.49154895e-09 

5.23979146e-07 

64 

5.43330535e-ll 

1.15360543e-08 

128 

8.48088179e-13 

2.46158319e-10 

256 

1.32479989e-14 

5.14359846e-12 

512 

2.06986786e-16 

1.05908706e-13 

1024 

3.2341 1700e- 18 

2.15760798e-15 

Table  1.  Screen  size  (N)  dependent  a-coefficient  scaling  factors. 


via  the  inverse  Fourier  transform:  <f>P(fc) — - — ><fip(r),  in  the  continuous  case  -  and,  in 

the  analogous  discrete  case:^p(m) — - — ></>P{p)-   This  will  be  accomplished  in  the 
discrete  two-dimensional  case,  according  to  the  following  equation:  (Fried,  Sep  97,  p.  4) 


<f>p{p)  =  N~l  ]T  ^(m)exp 


OT.7J=1 


ATtim  — 


\ 


N 


(15) 


where  m  =  (m,  n) .  Once  in  this  spatial  domain  form,  <j>p  (p)  will  consist  of  a  complex 

random  array,  the  real  and  imaginary  parts  of  which  each  will  correspond  to  a  random 

phase  pattern.  This  phase  pattern  is  suitable  for  high  spatial  frequencies  but  has  a 

deficiency  in  the  lowest  spatial  frequency  components.  Dr.  Fried  states, 

This  is  manifested  by  a  structure  function  for  this  pattern  which  has  the 
appropriate  five-thirds  power  law  dependence  for  separations  much 
smaller  than  the  real-world  physical  size  of  the  screen,  L,  but  deviates 
very  significantly  from  that  power  law  behavior  for  separations  as  large  as 
L/10,  and  even  decreases  in  value  as  the  separation  increases  beyond  L/2 
(Fried,  Sep  97,  pp.4-5). 

In  order  to   compensate  for  this  deficiency  in  the   lower  spatial   frequency 

components,  Fried  has  developed  a  suitably  chosen  random  tilt  term  that  can  be  added  to 
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(f>P(p).  The  results  with  this  added  term  prove  sufficient,  and  the  resultant  simulated 
phase  screen  will  be  accurate,  i.e.  its  structure  function  will  follow  a  five-thirds  power 
law  for  all  separations.  The  extra  tilt-term  can  be  seen  developed  and  validated  in  the 
reference  source  (Fried,  1996).  However,  in  this  treatise,  this  additional  tilt  term  is  not 
necessary  since  the  adaptive  optics  will  have  no  problem  canceling  that  part  of  the  phase 
variation.  This  particular  exercise  will  be  limited  to  laser  apertures  under  the  correction 
of  adaptive  optics,  and  as  such,  the  adaptive  optics  will  introduce  their  own  transfer 
function  -  eliminating  most  the  wavefront  distortion,  to  be  discussed  later.  The  pertinent 
point  here,  is  that  the  adaptive  optics  can  correct  the  tilt  terms  in  the  distorted  wavefront 
very  effectively,  thereby  eliminating  the  need  to  include  them  in  the  first  place  (Fried. 
1998). 

C.         ADAPTIVE  OPTICS  COMPENSATION 

As  mentioned  previously,  the  adaptive  optics  component  will  introduce  its  own 
Transfer  Function  (TF),  and  as  such,  will  require  treatment  in  this  simulation.  To  be 
precise,  the  adaptive  optics  will  correct  all  phase  fluctuations  to  the  extent  that  the  size  of 
the  sub-aperture  will  allow.  For  example,  if  it  were  possible  to  build  an  adaptive  optics 
system  with  infinitely  small  sub-apertures,  then  the  wavefront  distortion  could  be 
corrected  perfectly.  Allowing  for  finite  size  sub-apertures  there  will  be  residual  errors 
after  adaptive  optics  correction;  we  shall  consider  adaptive  optics  systems  with  finite  size 
sub-apertures,  thus  introducing  the  deviation  from  perfection  whose  statistics  are  studied 
here.  In  addition,  we  will  ignore  latency  effects  caused  by  the  time  delay  required  by  the 
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adaptive  optical  system.  The  imperfection  of  the  compensation  provided  by  the  adaptive 
optics  is  represented  by  what  is  called  the  residual  phase  error,  by  the  deviation  of  the 
corrected  wavefront  from  perfection;  this  is  called  fitting  error. 

In  this  simulation,  the  residual  fitting  error  will  be  developed  from  the  phase 
screen  through  the  use  of  its  transfer  function  mentioned  earlier.  As  developed  in 
Greenwood,  1979,  the  adaptive  optics  corrected  residual  phase  error  can  be  represented 
by  the  following  equation;  returning  to  a  familiar  continuous  treatment,  we  write 

0ao (f)  =  Jexp(2;r ifc-r)  TA0 (k)  0p (k)  die  ,  (1 6a) 


and 


$ao(k)=  Jexp(-2;rz*:-r)  (f>A0{r)dr 


(16b) 


=  Tao{k)(J)p{k) 
where  Ta0(k)  is  the  adaptive  optics  fitting  error  transfer  function  written  in  frequency 

terms.  <pp  (k)  is  the  Fourier  transfer  of  <j)p (r ) ,  which  is  the  phase  without  correction  from 
adaptive  optics,  as  previously  discussed.  As  shown  in  Greenwood,  1979,  the  value  of 
TA0 (ic)  is  given  by  the  following  equation: 


r*>(*K  l- 


2  J,  \7t\  K  \d) 

n\  K\d 


4J2(7r\ic\d) 
n\  K\d 


(17) 


where  d  denotes  the  effective  size  of  the  sub-aperture  components  in  the  adaptive  optics. 
As  can  be  seen  from  Figure  (4),  the  TA0  (k)  transfer  function  has  the  effect  of  attenuating 

the  lower  frequency  terms  and  passing  the  high  frequency  terms,  analogous  to  a  high-pass 
filter  in  electronics  (Greenwood,  1 979). 
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Figure  4.  Greenwood  Transfer  Function  Parallels  a  High-Pass  Filter 


STREHL  COMPUTATION 


The  computation  of  the  Strehl  ratio  turns  out  to  be  the  most  time-intensive  portion 
of  the  simulation.  However,  through  the  use  of  a  fast  Fourier  transform  (FFT)  and  clever 
assumptions,  one  statistically  random  phase  screen  realization  can  generate  many 
simulation  results.  In  the  case  of  one  phase  screen  represented  by  equation  (15)  and 
simulating  the  random  perturbations  a  light  ray  has  encountered  through  the  turbulent 
medium,  the  Strehl  ratio  can  be  computed  by  comparing  the  intensity  calculation  with  the 
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phase  screen  to  that  intensity  calculation  of  a  reference  intensity  with  zero  phase 
perturbations  (diffraction-limited). 

For  an  aperture  whose  center  is  at  r  ,  the  signal  power,  P(r) ,  comes  directly  from 
squaring  the  signal  strength,  S(r).  The  signal  strength  calculation  is  most  conveniently 
expressed  in  the  following  manner,  which  facilitates  future  development: 


r  -r  x 


S{r)=\A\—-\^V{i<f>p(r'))dF'i  (18) 

where   <f>P\j)   represents  the  previously  discussed  random  phase  distortion  function 
evaluated  at  f ' .  The  aperture  function,  A  (x) ,  is  defined  as  follows: 

A(x)=lh        '/|X"'"i      •  09) 

I  0,         if  otherwise 


Thus  A 


fr'-r\ 

represents  an  aperture  of  diameter  D  on  the  r  —  plane,  with  its  center  at 


D 

r'=r.  The  representation  provided  by  equation  (18)  is  beneficial  two-fold:  a)  the 
reference  intensity  is  easily  calculated  up-front  by  replacing  the  <f)P{r')  with  <f>P(r')  — 
zero,  and  b)  the  calculation  of  the  signal  strength  with  phase  distortion  and  adaptive 
optics  compensation  can  be  completed  easily  by  replacing  the  <f>P(r')  with 
(j)P(f')  =>  <fiAO(r')  from  equation  (16). 

It  is  here  that  Fourier  transforms  can  be  used  to  provide  a  tremendous  increase  in 
productivity.  As  described  in  the  previous  section  and  also  as  represented  in  equation 
(18),  it  is  assumed  that  the  f  vector  is  positioned  at  the  center  of  the  array  (although  not 
explicitly  stated  as  such)  and  that  the  r'  vector  moves  with  the  integration  along  all  of  the 
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points  over  the  aperture.   In  fact,  the  r  vector  does  not  need  be  positioned  at  the  center, 
and  as  Fourier  transforms  will  allow,  the  r  vector  can  be  positioned  at  every  point  and 
continue  with  the  integrations  through  the  virtual  use  of  the  two-dimensional  convolution. 
Following  Fried,  1997,  the  phase  exponential  will  be  defined: 

£/(r)  =  exp{i£,0(r)},  (20) 

so  that  equation  (18)  becomes: 


r  -r  ^ 


S(F)  =  JA\—-\U(r')dr'.  (21) 

Written  this  way  it  is  manifest  that  S(r)  is  equal  to  the  convolution  of  A  with  U.  A 
convolution  indicated  by  the  usual  (*)  designation,  of  the  two  elements  in  the  integrand  of 
equation  (21)  would  appear  as  follows: 

S(F)  =  A\'-^-)*U(F'),  (22) 

where  the  S(r)  result  would  now  be  an  Nx  N  matrix  of  the  same  dimensions  as  the  U(F') 
matrix.  Following  the  treatment  shown  by  Fried,  and  independently  applying  convolution 
manipulation,  equation  (22)  is  transformed  again  using  Fourier  transforms: 

S(r)=  \exp(2xifc-F)A{K)U(fc)dic,  (23) 

where  A  (k)  and  U  (re)  are  defined  by  the  following  equations: 

A(k)=  Jexp(-2;r/ic-r)  A{r)  dr.  (24a) 

U(k)=  jexp(-2xiK-f)U(r)df.  (24b) 
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One  other  point  of  note  concerns  the  value  of  the  Fourier  transformed  elements, 
and  in  particular  the  value  of  A(ic).  As  opposed  to  the  case  of  0(rc),  which  includes 
calculating  a  two-dimensional  FFT  of  the  randomly  generated  U(r),  the  value  of  A(k)  is 
a  closed-form  result  depending  on  the  geometry  of  the  aperture.  In  this  case,  a  uniform 
circular  aperture  is  used,  and  the  corresponding  result  can  be  expressed  as  follows: 


A(k)  =  -ttD2 
4 


2  J]  (tu\  K  |  D 
mK\D 


(25) 


This  expression  allows  a  more  refined  representation  in  the  computer  code,  avoiding  the 
computer  rounding  inaccuracies  that  arise  during  the  calculation  of  the  discrete  FFT.  A 
smooth  illustration  of  the  Fourier  transform  of  a  circular  aperture  is  shown  in  Figure  (5). 

E.         INCORPORATION  INTO  COMPUTER  CODE 

All  of  the  previous  sections  contain  the  mathematical  background  for  efficiently 
generating  many  simulations  of  Strehl  ratio  calculations  for  a  given  set  of  initial 
conditions.  As  will  be  detailed  later,  the  different  exponential  degrees  to  which  'many' 
can  cover  will  depend  on  the  speed  of  the  computer  being  used  as  well  as  how  many 
computers  (or  processors)  can  effectively  be  utilized  in  parallel.  This  section  will  present 
the  applicable  formulations  to  take  into  account  discrete  implementation  as  opposed  to 
the  continuous  treatment  previously  discussed. 

The  last  discrete  formulation  presented  was  Equation  (15)  representing  the  (j)P{p) 
phase  perturbation  term  in  frequency  space.  The  next  two  equations  present  the  discrete 
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Figure  5.  Fourier  Transform  of  a  Circular  Aperture  (Spatial  Frequency  Domain) 


version  of  the  adaptive  optics  compensation,  Ta0(k)  from  Equation  (16): 


x„  =jzdJi  k1  +  kI  , 

m  \         m  n     " 


and 


(26) 


TA0(m)=l\- 


2J,{xm) 


-i  2 


4«/2feJ" 


(27) 


Incorporating  Equation  (27)  into  the  resultant  phase  perturbation  and  performing  the 
inverse  discrete  Fourier  transform  result  in  the  following  two  equations: 


f 


\ 


1       ■  —         f 
-L7llKm  ■  — 

m    N 


\ 


and 


(28) 
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7V 

0Ao(P)=    Z    ^o(«)eXP 


ffi,n=l 


9      •-       ^ 

2niK_  •  — 


\ 


N 


(29) 


The  culminating  incorporation  into  the  discrete  adaptive  optics  compensating  phase 
screen,  UAO  (p) ,  is  given  by  the  anti-climactic  expression: 


UAo(p)=QXp{i<f>Ao(p)l 


(30) 


Its  discrete  Fourier  transformation  follows: 


UA0{m)=  ]T  UAO(p)exp 


'    0     -      P^ 

-L7ZIK„  ■  — 

m    Nj 


(31) 


In  order  to  discretely  represent  the  Fourier  transform  of  the  aperture,  a  conversion 
of  Equation  (25)  yields  the  following  two  discrete  representations: 

Xm=xD,]K>+KZ,      and  (32) 


1 

A(m)  =  -7uD' 
4 


2  J,  (AT.) 


X. 


(33) 


Possessing  discrete  formulations  for  both  UA0(m)  and  A(m),  the  discrete  conversion  of 


Equation  (23)  follows  easily: 

S(p)  =  N~'2ft  UA0(m)A(m)exp 


m,n=\ 


N  j 


(34) 


After  completion  of  Equation  (34)  in  the  computer  code,  we  are  effectively  left 
with  a  statistically  random  N  x  N  matrix  of  signal  strength  values,  all  for  the  same 
aperture  diameter  D  and  the  same  r0  value.  The  signal  strength  values  can  then  be 
squared  and  divided  by  the  reference  intensity  value,  yielding  an  N  x  N  matrix  of  Strehl 
ratio   values.   Figure   (6)   is  a  block  diagram  representation  of  the   functional   flow 
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throughout  the  computer  code.  It  starts  from  the  very  beginning  of  variable  initialization 
and  reference  intensity  computation  through  the  finishing  -  which  is  only  truly  finished 
when  the  desired  number  of  Monte  Carlo  realizations  has  been  reached.  Numbers,  for 
instance  through  the  completion  of  this  thesis,  range  from  two  loops  in  the  process  of 
testing  code  accuracy  to  83,000  loops  in  the  preparation  often  billion  realizations  for  one 
particular  r0  value.  A  later  version  of  the  computer  code,  written  in  Matlab  and  named 
'Sin.m'*,  is  provided  in  Appendix  (A).  As  such,  this  version  includes  extra  code  to 
account  for  intensity  scintillation  effects  to  be  discussed  later  in  chapter  IV. 

F.         RESULTS 

1.  Preparation 

Over  the  two-week  span  of  12-26  March  1998,  results  were  gathered  from 
implementation  of  the  aforementioned  computer  code,  which  at  the  time,  was  named 
'turro.m/  The  following  parameters  will  be  analyzed  as  a  result  of  those  computer  runs, 
although  all  may  be  changed  to  run  any  particularly  interesting  case.  An  aperture  diameter 
size  of  1.2  meters  will  represent  the  telescope,  with  adaptive  optics  sub-aperture 
components  having  an  effective  diameter  size  of  1 5  centimeters.  The  resultant  D/d  ratio 
in  this  case  becomes  8.0.  The  simulated  size  of  the  phase  perturbation  screen  will  be  3.0 
meters,  thereby  fully  covering  the  simulated  aperture.  Finally,  the  range  of  eleven  r0 
values  will  be  applied,  from  12.6  to  14.1  centimeters  in  increments  of  0.015  centimeters. 


'Sin'  is  short  for  'scintillation'  in  this  particular  case,  not  to  be  confused  with  any  religious  meaning. 
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Figure  6.  Block  Diagram  of  Functional  Flow  in  the  Computer  Code 

The  resultant  d/r0  ratios  in  this  case  become  1.190  down  to  1.064.   These  parameters  are 
summarized  in  Table  (2). 
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Parameter 

Initialized  Value 

Aperture  Diameter 

1.2  (m) 

Sub-Aperture  Diameter 

0.15  (m) 

D/d  Ratio 

8.0 

Screen  Diameter 

3.0  (m) 

Values  ofr0  (meters) 

0.1410,  0.1395,  0.1380,  0.1365,  0.1350. 
0.1305,0.1290,0.1275,0.1260 

0.1335, 

0.1320, 

Resultant  d/r0  Ratios 

1.0638,  1.0753,  1.0870,  1.0989,  1.1111, 
1.1494,1.1628,  1.1765,  1.1905 

1.1236, 

1.1364, 

Values  of  Scintillation 

Zero 

Screen  Size(N) 

256x256 

Table  2.  Initialized  Parameters  for  Results  of 'Turro'  Computer  Runs,  12-26  March  1998. 

As  mentioned  earlier,  the  design  of  the  computer  code  allows  many  simulations  to 
be  run  on  different  computers,  effectively  parallel  processing  without  actually  linking  the 
processors.  In  fact,  the  storing  of  the  values  in  a  sparse  matrix  (histogram  fashion), 
allows  the  programmer  to  gather  the  results  from  each  separate  computer  that  was  utilized 
and  add  all  of  the  results  together  into  one  combined  sparse  matrix.  Figure  (7)  below 
displays  the  results  that  were  generated  from  the  above  initialized  values.  The 
observation  of  "830  million  per  line"  indicates  the  number  of  Monte  Carlo  simulation 
results  that  built  each  apparently  Gaussian  shaped  bell  curve.  Each  run  through  the  main 
loop  of  the  program  yields  2  x  256  or  131,072  realization  results.  One  quick  calculation 
reveals  that  6336  runs  through  the  main  loop  were  needed  for  each  r0  value,  yielding  a 
total  of  69696  runs  for  all  eleven  values  of  turbulence  severity. 

To  appreciate  the  magnitude  of  this  number,  consider  the  following  actual 
numbers  obtained  during  the  generation  of  the  lines  on  Figure  (7).  One  run  through  the 
main  loop  of  the  computer  code,  executed  in  Matlab  on  a  200  MHz  Intel®  Pentium®  Pro 
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Figure  7.  12-26  March  1998  'Turro'  Computer  Code  Results 

(circa  1997)  personal  computer  (PC),  requires  approximately  one  minute.  The  required 
run  time  would  result  to  484  days  which  is  clearly  well  over  a  year  for  completion. 
Utilizing  access  to  a  multitude  of  Sun®  Sparc®- workstations  (circa  1990)  across  the 
Naval  Postgraduate  School  campus  proved  to  be  the  saving  grace.  Although  the  time  to 
complete  the  main  loop  on  any  individual  Sun  was  significantly  slower  than  the  PC,  the 
greater  number  was  invaluable. 

A  familiarity  with  the  Unix®  operating  system  is  essential  to  efficiently  gather 
data  from  multiple  runs  across  the  campus  network.  One  needs  first  to  check  the 
prospective  machine  to  determine  the  processes  already  running  and  then  start  the  Matlab 
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programs  running  in  the  background  before  moving  to  the  next  machine.  The  'Tsh" 
command  is  needed  to  log-on  quickly,  complete  business  on  the  machine,  and  then  exit  to 
"rsh"  on  another.  At  the  time  of  execution  in  March  of  1998,  there  were  seventeen 
Sparc®-stations  in  Ingersol  Hall,  thirteen  in  Root  Hall,  and  eight  in  Spanagel  Hall  for  a 
total  of  38  prospective  machines.  Only  28  of  the  machines  turned  out  to  be  'available" 
after  checking  to  ensure  that  nobody  else  was  running  their  own  programs,  taking  up 
precious  processor  usage  time.  Over  the  course  of  nearly  a  week  on  28  machines,  the  830 
million  simulations  were  realized. 

2.  Examination  of  Results 

The  data  displayed  in  Figure  (7)  was  gathered  in  such  a  way  as  to  allow 
examination  in  the  following  manner.  In  order  to  analyze  low-probability  events,  the  data 
must  be  converted  from  a  histogram  showing  a  bell  curve  into  a  corresponding 
probability  distribution  graph.  The  histogram  is  essentially  a  probability  density-  curve  for 
each  curve.  Computing  the  cumulative  probability  is  necessary  in  order  to  generate  the 
next  graph,  displayed  as  Figure  (8).  A  curious  circumstance  of  working  in  such  extreme 
low  probabilities  is  the  number  of  'nines'  (high),  or  'zeros'  (low).  A  probability  of 
0.0000001  =  10"7  is  displayed  in  the  graph's  y-axis  as  (-7),  short  for  1.0E-7,  and 
conversely  a  probability  of  0.9999999  =  1-10"7  is  displayed  as  (7),  short  for  1.0E-7  or  7- 
nines  (high). 

The  next  order  of  business  is  that  of  prediction  —  or  at  the  very  least,  is  there 
some  natural  mathematical  law  that  is  discernable  from  the  data?  Referring  to  Figure  (8) 
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Figure  8.  Strehl  vs.  Probability  with  Crossing  Line 

again,  there  is  a  vertical  broken  line  drawn  at  Strehl  =  0.30.  One  can  visualize  an  infinite 
number  of  parallel  curves  above  and  below  the  eleven  curves  shown,  and  the  vertical  line 
graphically  cuts  through  all  eleven  and  extends  outward  into  the  infinite  number  of  other 
such  curves  not  shown,  nor  even  generated.  Changing  the  axes  of  the  graph  and  plotting 
differently,  Figure  (9)  below  shows  five  such  vertical  lines  plotted  as  d/r0  versus 
probability.  The  lines  appear  to  be  very  straight.  Nevertheless,  exactly  how  straight  are 
they? 

To  check  the  straightness  of  the  lines,  it  is  necessary  to  generate  a  separate  set  of 
data  at  less  severe  turbulence  values,  or  higher  r0  values.  Take  the  higher  Strehl  =  0.370 
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Figure  9.  Turbulence  Severity  vs.  Probability  for  selected  constant  Strehl  lines 

line  for  example  from  Figure  (9).  One  can  "fit"  an  exact  line  to  the  data  points  that  are 
plotted  as  an  apparent  straight  line.  That  same  Strehl  =  0.370  fitted  line  can  then  be 
extrapolated  further  to  the  bottom-left  of  the  graph.  Figure  (10)  displays  that  detailed 
extrapolation  for  the  top  three  lines  from  the  previous  Figure  (9).  Figure  (10)  also 
displays  the  results,  shown  as  small  circles  according  to  the  legend,  of  Strehl  probability 
data  generated  separately  according  to  the  initialized  variables  shown  in  Table  (3).  The 
initialization  values  shown  in  Table  (3)  were  chosen  to  provide  "checker"  data  points.  As 
plainly  seen  from  Figure  (10),  the  "check"  data,  labeled  as  "test"  data,  does  not  fall  on  the 
extrapolated  line  from  the  "fit-to"  data.  At  the  time,  this  was  a  very  discouraging  result. 
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Parameter 

Initialized  Value 

Aperture  Diameter 

1.2  (m) 

Sub-Aperture  Diameter 

0.15  (m) 

D/d  Ratio 

8.0 

Screen  Diameter 

3.0  (m) 

Values  ofr0  (meters) 

0.1590,  0.1560,  0.1530,  0.1500,  0.1470,  0.1440 

Resultant  d/r0  Ratios 

0.9434,0.9615,0.9804,  1.0000,  1.0204,  1.0417 

Values  of  Scintillation 

Zero  compensation 

Screen  Size(N) 

256  x  256 

Table  3.  Initialized  Parameters  for  "Test"  Data  of 'Turro'  Computer  Runs,  12-26  March  1998. 

Evaluation  that  is  more  extensive  led  to  the  following  question.  If  not  straight  on  a  linear 
scale,  how  straight  is  the  line  on  a  log  scale? 
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Keeping  the  dependent  axis  the  same,  but  changing  the  independent  axis  to  a  log 
scale  proved  to  be  the  answer.  Figure  (11)  shows  the  results  from  such  a  change.  The 
fitted  line  was  fitted  to  the  logarithm,  base  10,  of  the  "fit-to''  data.  The  same 
extrapolation  was  performed  and  in  the  same  manner,  the  small  circles  represent  where 
the  "checker"  data  resulted.  As  can  also  be  seen  plainly  from  Figure  (11),  the  data  fall  on 
a  very  straight  line,  indeed,  albeit  a  straight  logarithmic  line. 
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IV.      STREHL  DEGRADATION  DUE  TO  PHASE  AND 
SCINTILLATION  EFFECTS 


A.         OVERVIEW 

It  is  a  fact  that  adaptive  optic  systems  are  indeed  proving  very  effective  by  just 
compensating  for  phase  effect  distortions  along  the  path  of  light  propagation.  However,  it 
is  also  considered  interesting  and  possibly  extremely  helpful  to  model  the  effects  of 
intensity  scintillation  and  combine  the  model  with  that  already  discussed  in  an  attempt  to 
evaluate  the  same  low  probability  Strehl  events.  As  mentioned  previously,  the  effects  of 
intensity  scintillation  can  be  simulated  separately  and  in  the  same  manner  as  the  phase 
effects,  that  of  finding  an  expression  to  represent  the  power  spectral  density.  Once  that 
expression  is  determined  and  then  validated,  it  can  be  included  into  the  'turro.m'  program 
to  again  generate  many  simulations  for  a  given  r0  value  and  a  given  value  of  intensity 
scintillation,  af  . 

For  the  following  development,  it  is  assumed  that  the  particular  optical  system  in 
question  has  the  capability  to  perfectly  correct  phase-only  perturbations  with  its  adaptive 
optics  components.  We  will  continue  to  use  the  principle  of  reciprocity  to  analyze  a  laser 
beam  propagating  outward  and  upward  though  the  atmosphere  towards  a  target  in  space. 
And  as  before,  we  will  only  be  concerned  with  an  evaluation  of  the  statistics  of  this 
scintillation  effect.  As  such,  and  following  Fried  (from  Fried,  1997,  p. 3),  the  statistical 
notation,  i(F),  called  the  turbulence  induced  random  log-amplitude  at  position  r,  will 


represent  all  of  the  effects  of  turbulence  that  will  be  of  concern  to  us  in  this  thesis.   The 
next  section  will  address  the  details  of  its  development. 

B.         SCINTILLATION  SCREEN  INCORPORATION 

The  purpose  of  this  section  is  to  develop  a  similar  random  screen  to  the  phase 
perturbation  screen  from  the  previous  chapter,  with  this  particular  screen,  <f>( ,  having  the 

statistical  characteristics  of  atmospheric  turbulence  induced  intensity  scintillation.   More 
explicitly,  we  wish  to  able  to  express  an  optical  field,  U(r),  in  the  following  expression: 

tf(r)  =  V2expfo(r)).  (35) 

This  will  prove  desirable  as  we  can  again  use  equation  (21),  developed  earlier  in  chapter 
three  and  reproduced  here  for  the  benefit  of  the  reader: 

S{¥)=\A["—^\u{V')dF'.  (21) 

1.  Log- Amplitude  Fluctuation  Statistics 

The  random  log-amplitude,  £ (r) ,  as  developed  in  the  turbulence  literature  from 
Kolmogorov  to  the  present,  may  be  considered  to  be  a  random  function  characterized  by  a 
mean  value  1  and  a  covariance  function  C(  (?) .  These  two  functions  are  defined  by  the 
following  equations:  (Fried,  1997,  p. 5) 

1  =  (£(F)},  (36) 

and       Ce(r)=([t(r'+±r)-l][£(F'-±r)-!])  .  (37) 
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A  very  instrumental  parameter  follows,  called  the  log-amplitude  variance  and 
denoted  by  ae2.  It  corresponds  to  the  r  =0  value  of  C(  (r) ,  which  is  to  write  explicitly: 

cr,2=C,(0).  (38) 

The  importance  of  this  parameter  stems  from  the  following  relationship  that  can  be 
shown  from  conservation  of  energy  principles: 

!  =  -o).  (39) 

This  states  simply  that  given  a  value  of  g(2,  we  empirically  know  the  mean  value  of  the 
random  log-amplitude  function  (Fried,  1997,  p. 5). 

The  following  equation,  (40),  follows  from  more  detailed  development  presented 
in  Fried,  1997.  The  pertinent  points  to  be  taken  from  this  development  include  the  fact 
that  the  Ce  (r)  function  depends  on  the  path  and  the  particular  wavelength.  The 
propagation  path  in  this  thesis  runs  from  s-0  on  the  ground  to  s  =>  oo  at  the  space  end 
of  the  path.  The  notation  k  again  refers  to  the  optical  wavenumber  related  to  the 
wavelength  by  the  expression  k=2nl X  .  Presented  without  derivation,  we  write:  (Fried, 
1997,  p.6) 

^=Bk1/6    J  Cl(s)  s5''6  ds  ,  (40) 


PATH 


where  the  value  of  B  is  given  as  B*  0.56315761 .    Hence,  we  are  left  with  a  means  to 
measure,  or  in  this  case,  generate  a  given  value  of  <5(  . 

The  other  important  parameter,  the  covariance  function  C(  (F) ,  is  needed  in  order 

to  generate  the  power  spectral  density,  dV   Once  given  a  power  spectral  density,  we  can 


impart  that  statistic  onto  the  random  generated  values  representing  the  spatial  frequencies 
-  exactly  as  done  in  the  previous  chapter  when  producing  the  phase  perturbation  screens. 
The  important  results  will  again  be  presented,  allowing  the  interested  reader  to  review  the 
source  document  for  details.  After  appreciable  development  in  Fried  (1997,  pp. 6-7),  we 


write: 


Ce(r)  =  jexp(2/r/A:  -r)  <&((k)  dK  , 


(41) 


where  it  is  noted  that  the  two  functions  are  Fourier  transform  pairs.  The  power  spectral 


density,  O^,  is  given  by  the  equation: 


0,(*:)  =  A£2(2;rr8/3|^|-11/3    \  C2N(s) 


PATH 


f 


1-cos 


(2tt)2  k 


\ 


ds , 


(42) 


where  the  value  of  A  is  given  as  A*  0.65150030  and  the  path  is  the  same  as  above. 


2. 


Screen  Formulation 


Incorporation  into  scintillation  screen  computer  code  follows  very  closely  to  that 
presented  in  the  previous  chapter.  The  generation  of  the  array  of  random  values  is 
accomplished  in  the  same  manner,  by  defining  the  spatial  frequency  components  of  the 
array,  w(m).  These  elements  will  again  be  complex  and  take  on  random  values  with  zero 
mean  and  unity  variance.  For  the  sake  of  clarity,  the  two  components  of  the  m  vectors 
will  be  Km  and  Kn ,  with  the  values  before  randomization  initialized  as  follows: 


^rnandKn  m  ™/L, 


(43) 


where  m  will  be  defined: 
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^      f(w-l)  if       \<m<\N 

m  =  \  .  (44) 

\{m-N-\)       if      \N<m<N 

The  second  step  will  also  be  the  same  as  the  random  values  are  scaled  to  be 
proportional  to  the  magnitude  of  the  square  root  of  the  log-amplitude's  power  spectral 
density.  Corresponding  to  equation  (12)  in  the  earlier  development,  this  can  be  expressed 
in  the  spatial  frequency  domain  as  follows: 

l{m)  =  w(m){a^®e{m)  )    ,  (45) 

where  the  a  scaling  factor  is  defined  as  a  =  N2/L  in  order  to  account  for  the  physical 
size  of  the  screen.  This  will  be  apparent  in  the  next  step,  where  the  result  of  the  inverse 
Fourier  transform,  as  realized  by  an  inverse  FFT,  will  now  produce  a  suitable  magnitude 
result.  (Fried,  1997,  pp.  10-11) 

The  third  step  corresponds  to  the  recently  developed  equation  (40),  replacing  the 
integration  with  a  summation: 

N   ~  (  m\ 

£(p)  =  N-2  Y,  ^)exp   2ni—     .  (46) 


m,n=\ 


K 


For  the  sake  of  clarity,  the  £  (p)  result  is  a  two-dimensional  complex  array  of  random 
values  denoting  spatial  coordinates,  remembering  prior  to  the  final  step  that  the  mean 
value  here  is  still  zero. 

All  that  is  left  to  do  in  the  fourth  and  final  step  is  to  apply  the  mean  value  as 
defined  in  equations  (36)  and  (39).  This  is  realized  by  computing  equation  (40)  in 
discrete  increments  (i.e.  estimating  the  integration  by  selecting  a  finite  number  of 
intervals  in  the  path).    The  mean  value  adjustment  is  made,  yielding  the  screen(s)  and 
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written  in  the  following  two  relationships: 

Mp)  i  =*{'(p)}-0-?,  (47) 

and       MP)2=*{t(.P)}-<rl    •  (48) 

It  is  at  this  point  we  have  a  statistically  sound  representation  of  the  log-amplitude 
function,  and  the  interested  readers  are  referred  to  the  source  document  (Fried,  1 997,  pp. 
11-15)  for  verification  of  the  accuracy  with  this  method.  These  screens  represented  by 
equations  (47)  and  (48)  are  now  ready  to  implement  in  the  computer  code  as  an  extension 
to  equation  (30),  written  as  follows: 

UM+l(p)  =  &v{iJAo(p)Yexp{MP)}  ■  (49) 

There  is  no  deviation  in  the  remainder  of  the  'turro'  computer  from  this  point  on. 

A  word  about  the  real-world  practical  significance  of  the  turbulence  induced  log- 
amplitude  value  is  relevant  here.  The  unit  of  measurement  for  oe  is  the  square  Napier, 
which  is  denoted  as  Np  .  Actual  turbulent  airflow  measurements  range  from  an  ideal  of 
zero  through  a  typical  measurement  of  oe  =  0.05  Np  ,  up  to  a  severe  turbulent  condition 

9  9  9 

of  Ge  =  0.15  Np  .  It  is  with  this  g(  value  that  the  turbulence  severity  will  be  represented, 
just  as  the  r0  value  represented  the  turbulence  severity  in  the  phase  perturbation  of  the 
beam. 

C.        RESULTS 

Over  the  two-week  span  from  21  October-  4  November  1998,  results  were 
gathered  from  implementation  of  the  combined  computer  code,  which  at  this  time  is 
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named  'sin.m.'  The  following  parameters  will  be  analyzed  as  a  result  of  these  computer 
runs,  although  just  as  earlier  mentioned,  all  initialized  parameters  may  be  changed  to  run 
any  particularly  interesting  case.  An  aperture  diameter  size  of  1 . 1  meters  will  represent 
the  telescope,  with  adaptive  optics  sub-aperture  components  having  an  effective  diameter 
size  of  10  centimeters.  The  resultant  D/d  ratio  in  this  case  becomes  1 1.0.  The  simulated 
size  of  the  phase  perturbation  screen  will  be  3.0  meters,  which  again  will  fully  cover  the 
simulated  aperture. 

With  additional  code  to  account  for  the  intensity  scintillation,  the  average  time  to 
complete  one  main  loop  of  the  code  is  increased  to  approximately  one  minute  thirty 
seconds  on  the  same  200  MHz  PC.  Consequently,  the  range  of  r0  values  will  be  reduced 
to  three,  namely  10.2,  11.1,  and  12.0  centimeters.  The  resultant  d/r0  ratios  in  this  case 
become  0.980  down  to  0.833 .  Finally,  each  of  the  r0  values  will  be  coupled  with  five 
intensity  scintillation  values  ranging  from  a  very  severe  o^2  =  0.1 1  Np2  all  the  way  down 
to  zero  scintillation.  These  parameters  are  summarized  as  well  below  in  Table  (4). 


Parameter 

Initialized  Value 

Aperture  Diameter 

1.1  (m) 

Sub-Aperture  Diameter 

0.10  (m) 

D/d Ratio 

11.0 

Screen  Diameter 

3.0  (m) 

Wavelength  (X) 

0.840  (microns) 

Values  ofr0  (in  meters) 

0.1200,0.1110,0.1020 

Resultant  aVr0  Ratios 

0.8333,  0.9009,  0.9804 

Values  of  Scintillation  (in  Np  ) 

0.11,0.07,0.03,0.01,0.00 

Screen  Size(N) 

256x256 

Table  4.  Initialized  Parameters  for  Results  of 'Sin*  Computer  Runs,  20-31  October  1998. 
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The  long  time  delay  between  the  two  sets  of  runs  (seven  months)  proved  to  be 
beneficial  in  one  respect:  it  allowed  new  workstation  computers  to  be  acquired,  greatly 
improving  processor  speeds  and  memory  availability.  The  following  histogram  results, 
displayed  in  Figure  (12)  below,  represent  two  billion  simulated  realizations  of  the 
combined  computer  code.  The  results  are  also  displayed  with  a  legend,  clearly  illustrating 
the  gradual  Strehl  ratio  decline  for  each  of  the  three  r0- values  as  the  intensity  scintillation 
gets  increasingly  worse. 

The  long  time  delay  between  the  two  runs  also  proved  to  be  detrimental  in  another 
respect:  the  results  were  actually  gathered  after  the  author's  graduation  date  and  thereby 
represented  a  successful  stopping  point  in  the  life  of  the  thesis.  Time  is  not  available  in 


x  10 


Strehl  Distribution 


Strehl  bin  width  =  0.0001 
2.06  x  109  per  curve 
D/d  =  11.0 


d/r0=  0.833 
d/r0=  0.901 
d/r0=  0.980 


Figure  12.  20-31  October  1998  'Sin'  Computer  Code  Results 
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this  thesis  for  detailed  examination  of  the  results,  as  we  were  able  to  do  in  the  previous 
chapter.  Only  cursory  remarks  will  be  made  during  the  following  presentation  of  the 
various  'snapshot'  figures,  with  the  intention  of  stimulating  thought  for  further 
examination  in  another  thesis,  possibly.  The  cursory  examination  will  be  presented  in 
two  segments.  The  first  segment  is  shown  as  "blow-up"  views  of  each  of  the  three 
separate  r0- value  groups,  followed  by  the  five  views  of  each  of  the  a,2-value  groups. 

Close  evaluation  of  Figure  (12)  reveals  one  very  evident  feature:  the  bell-shaped 
characteristic  of  each  curve  gets  wider  for  each  increase  in  a,2.  An  even  clearer 
manifestation  of  this  is  seen  in  Figures  (13),  (14)  and  (15)  where  it  can  be  seen  that  the 
slope  of  the  curves  decreases  as  af2  increases.  The  rotational  differences  in  the  slope  of 
the  probability  lines  towards  the  clockwise  direction  are  indicative  of  the  widening. 
There  is  also  an  appearance  that  the  lines  would  actually  cross  high  and  to  the  right,  if  in 
fact,  the  lines  were  extended  further  to  the  right. 

We  will  make  two  remarks  here  concerning  the  previous  observation.  The  first 
remark  is  a  confirmation  of  earlier  hypothesis  and  very  limited  actual  measurement.  The 
af2-value  has  a  greater  impact  on  the  variance  of  the  Strehl  ratio  measurement.  The 
decrease  in  slope  as  the  intensity  scintillation  increases,  so  indicated  in  the  three  previous 
figures,  is  in  direct  support  of  the  increase  in  variance.  Statistical  calculations  on  each  of 
the  bell  curves  support  an  increase  in  both  the  standard  deviation  and  the  variance  as  the 
intensity  scintillation  increases. 
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Figure  13.  Strehl  vs.  Probability  Lines  at  d/r0  =  0.833 
The  second  remark  will  require  extra  work  in  the  future  for  validation.  This 
concerns  the  overall  r0-value  after  both  r0  and  af  parameters  have  been  combined.  For 
example:  we  first  selected  an  r0-value  (i.e.,  r0  =  12  cm)  and  then  selected  a  scintillation 
value  (<j(  =0.05  Np  ),  without  readjusting  the  r0-value  for  the  worsening  that  will  occur 
as  a  result  of  the  addition  in  scintillation.  In  fact,  the  actual  r0-value  would  worsen  (to 
approximately  r0  =  10.5  cm)  because  the  added  intensity  scintillation  has  an  effect  on  the 
r0- value  measured  at  the  aperture.  An  approximate  solution  would  be  to  centroid  align  a 
set  of  intensity  variance  curves  leaving  d/r0  constant.  Unfortunately,  since  the  shapes  of 
the  histograms  change  slightly  for  different  d/r0,  this  is  not  exact.  Further  research  will  be 
needed  to  resolve  this  interaction. 
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Figure  14.  Strehl  vs.  Probability  Lines  at  d/r0  =  0.901 
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Figure  15.  Strehl  vs.  Probability  Lines  at  d/r0  =  0.980 
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It  is  the  author's  contention  that  as  the  intensity  scintillation  increases,  the  starting 
r0-value  would  need  to  be  weaker  (i.e.,  r0  =  13.5  cm),  so  that  after  applying  the  same 
scintillation  value  (<je  =  0.05  Np  ),  the  desired  r0- value  would  deteriorate  (to  r0  =  12  cm). 
This  'fix'  would  transpose  each  of  the  lines  to  the  right  in  proportion  to  the  increase  in 
oe  ,  forcing  the  intersection  to  move  down  to  the  center  position  of  0.5  probability  - 
effectively  having  all  five  of  the  lines  rotate  around  one  central  position. 

The  remaining  five  figures  (16-20)  present  the  same  information  organized  in  a 
different  manner.  In  examining  these  figures,  it  will  be  noted  that  the  Strehl  ratio  shifts  to 
the  left  as  each  figure  represents  a  constant  <je  -value  with  varying  phase-only  r0-values. 
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Figure  16.  Strehl  vs.  Probability  Lines  for  <5e   =  0.11  Np" 
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Figure  17.  Strehl  vs.  Probability  Lines  for  <Je   =  0.07  Np2 
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Figure  18.  Strehl  vs.  Probability  Lines  for  af  =  0.03  Np: 
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Figure  19.  Strehl  vs.  Probability  Lines  for  a/  =  0.01  Np2 
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Figure  20.  Strehl  vs.  Probability  Lines  for  a/  =  0.0  Np2 


50 


V.       CONCLUDING  REMARKS 

A.         LESSONS  LEARNED 

The  list  of  the  many  lessons  that  were  learned  will  be  shortened  to  include  only 
the  most  important  ones.  The  big  lessons  that  must  be  learned  very  well  include  the 
following: 

D    Become  intimately  familiar  with  the  Matlab  programming  language. 

D    Learn  the  basics  of  the  Unix  OS,  well  enough  to  movefiles  around  (using  ftp 
and  telnet),  and  execute  files  in  the  background. 

D    Get  access  to  the  fastest  workstation  processors  available  that  can  run 

Matlab.  At  the  time  of  this  printing,  on  the  NFS  grounds,  this  includes  the  SGI 
Octane®  and  the  Sun  Ultra  60  Sparcstation®. 

Smaller  lessons  are  important  as  well,  because  they  will  often  be  the  cause  of 
pounding  headaches.  First,  it  is  imperative  to  have  one  different  program  for  each 
separate  machine,  in  order  to  save  to  separate  data  files.  The  separate  data  files  will 
continue  to  grow  as  the  program  progresses,  and  when  finished,  the  data  files  can  be 
loaded  one-by-one  into  Matlab  and  recombined  by  adding  together  the  individual  sparse 
matrices.  Second,  the  random  number  generator  must  be  well  understood.  A  key 
statement  that  should  be  inserted  into  each  file  is  the  following: 
"randn('state',sum(100*clock))."  Otherwise,  Matlab  will  begin  each  individual  session 
with  the  exact  same  'state'  and  then  follow  with  the  same  random  numbers  every  time  — 
it  is  known  as  a  pseudo-random  number  generator. 
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During  the  course  of  this  thesis,  several  problems  were  encountered  while 
attempting  to  run  Matlab  scripts  in  the  background  on  different  workstations  around  the 
NPS  grounds.  First  of  all,  the  Matlab  code  must  be  a  script  to  run  in  the  background. 
Functions  will  not  process.  Good  information  can  be  found  on  this  topic  at  the  Matlab 
web  site.  Next,  it  is  highly  recommended  to  check  the  processes  running  on  a  particular 
machine  -  before  starting  a  Matlab  process  in  the  background.  The  'ps'  command  will 
perform  this  function,  with  the  'ps  -aux'  switches  working  well  on  SunOS  systems,  and 
'ps  -aef  '  working  well  on  Solaris  and  IRIX  systems.  Anything  less  than  99  percent 
processor  usage  time  is  detrimental  to  the  speed  of  processing  one's  program,  as  well  as 
increasing  the  chances  of  death  by  interference.  Working  over  long  weekends  and  holiday 
breaks  is  a  good  practice  to  follow. 

B.         POSSIBLE  FUTURE  RESEARCH 

This  would  best  be  stated  as  "where  do  we  go  from  here?"  There  are  several 
avenues  of  research  that  can  be  pursued.  Just  as  the  scintillation  was  added  to  the 
simulation  from  the  initial  effort,  more  effects  can  be  modeled  and  effectively  added  to 
the  simulation.  Optimization  of  the  computer  code  so  that  it  may  be  run  on  faster 
mainframe  computers  is  an  interesting  field  of  work  as  well.  This  section  will  address 
just  a  few  of  these  avenues,  and  rest  assured,  close  liaison  with  Dr.  Fried  is  strongly 
advised  (vital)  for  any  work  surrounding  this  computer  code.  This  is  especially  true  now 
that  it  has  been  changed  quite  a  bit  by  the  author. 
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Within  the  operational  field  of  adaptive  optics,  there  are  two  especially  promising 
topics  thought  to  be  "doable"  with  respect  to  realistic  simulation  on  computers.  The  first 
topic  is  called  branch  points.  Branch  points  are  manifested  in  actual  field  measurements 
from  adaptive  optics  systems,  specifically,  when  attempting  to  correct  phase  perturbations 
in  the  presence  of  strong  scintillation.  Very  simply,  adaptive  optic  corrections  are 
'measured'  by  sensing  the  phase  of  adjacent  light  rays  in  the  focal  plane,  and  then 
calculating  differences  between  the  measurements.  When  there  is  zero  intensity  in  one 
particular  spot  on  the  focal  plane,  there  is  a  2n  anomaly  in  the  phase  which  the  adaptive 
optics  system  cannot  cope  with  properly.  This  is  known  as  a  branch  point  from 
mathematics  theory  and  it  causes  a  malfunction  in  the  sensor.  Dr.  Fried  has  some 
thoughts  on  how  this  might  be  represented  in  this  type  of  computer  simulation. 

One  other  promising  topic  is  called  latency  within  the  adaptive  optics  field. 
Inherent  in  the  process  of  phase  corrections  is  the  fact  that  a  correction  is  inevitably  being 
applied  based  on  a  measurement  taken  a  short  time  (on  the  order  of  milliseconds)  prior.  It 
should  be  possible  to  develop  an  algorithm  capable  of  making  a  more  informed 
correction,  based  not  only  previous  data,  but  also  an  educated  (algorithmic)  prediction. 
Successful  efforts  in  this  area  would  reduce  feedback  problems  associated  with  the 
sensors  and  enable  higher  Strehl  ratio  capabilities. 

There  are  three  more  topics  that  present  challenging  difficulties,  but  not  much 
thought  has  been  devoted  to  solutions.  They  are  listed  briefly  here,  intended  only  to 
stimulate  thought: 
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□  Measurement  noise, 

□  Non  linear  signal  processing,  and 

D    Servo  effects  not  covered  by  a  time  delay. 

With  respect  to  computer  code  manipulation,  the  need  for  speed  in  processing 
time  never  diminishes.  Short  of  waiting  for  the  next  generation  faster  processor  and 
faster  memory,  the  only  other  way  is  to  "streamline"  the  code.  During  the  course  of  this 
thesis,  we  attempted  to  use  a  Matlab  compiler  with  the  intention  of  generating  an 
executable  file  in  C  or  C++.  Conversion  into  executable  code  would  enable  adaptation 
onto  much  faster  mainframe  computers,  such  as  SGI  Cray  machines,  for  example. 
Unfortunately,  there  were  two  important  functions  within  Matlab  ('sparse'  and  'spline') 
not  supported  by  the  Matlab  compiler  -  curious  as  that  was.  We  also  considered 
converting  the  entire  program  into  C++,  but  quickly  found  that  daunting  task  more  suited 
to  a  team  of  C++  programming  specialists.  Recently,  newer  versions  of  the  Matlab 
compiler  are  being  released  with  claims  of  support  for  previously  unsupported  functions, 
and  would  thereby  be  worth  trying. 

One  other  important  point  concerning  Matlab  that  may  be  important  when  porting 
to  multi-processor  machines  is  the  limitation  that  it  cannot  utilize  them  efficiently. 
Fortunately  for  this  application,  the  parallel  nature  of  the  code  not  only  allows  multiple 
runs  to  process  on  separate  machines,  but  multiple  runs  capable  on  the  same  multi- 
processor machine.  Each  of  the  processors  within  the  computer  will  take  dedicated  time 
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for  its  respective  program,  and  there  is  no  loss  in  processing  time  -  effectively  getting 
twice  the  output  in  the  same  amount  of  time  on  a  dual-processor  machine. 

C.         CONCLUSION 

In  this  thesis,  over  two  billion  computer-simulated  realizations  verified  a  log- 
normal  relationship  between  the  Strehl  ratio  and  its  probability  distribution  value  for  a 
given  d/r0  value,  and  similarly  a  log-normal  relationship  between  the  d/r0  value  and  its 
probability  distribution  value  for  a  given  Strehl  ratio  value.  Extremely  low  probabilities, 
down  to  1.0  x  10"  ,  correspond  to  the  number  of  realizations  collected,  with  stability  in 
the  lines  down  to  a  probability  of  5.5  x  10"6.  Most  recently,  over  the  course  of  a  week,  21 
computer  workstations  plus  another  5  dual  processor  workstations  (effectively  acting  as 
10)  were  running  31  separate  copies  of  the  computer  code  to  generate  15  curves. 

This  computer  code  is  available  to  run  simulations  for  any  given  adaptive  optics 
aperture  arrangement  in  any  given  atmospheric  turbulence  condition.  Through  the 
principle  of  reciprocity,  the  code  is  applicable  to  laser  beam  transmission  up  through  the 
atmosphere,  as  well  as  light  collection  in  an  astronomical  telescope.  It  can  be  scaled  to 
run  any  desired  number  of  realizations,  given  the  proper  computer  resources  and  realistic 
time  constraints.  The  histogram  results  are  conducive  to  probability  analysis  and  can 
therefore  be  used  in  conjunction  with  other  criteria  to  make  informed  decisions  about 
prospective  adaptive  optics  system  design 
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APPENDIX  A.    MATLAB  COMPUTER  CODE  'SIN.MP 


Sin.m 

'/. 

'/.    Scripti  called  Sin  (short  for  scinti  llation)  i  to  incorporate  Dr.  Fried's 

'/.    programs-   Both  of  these  function  programs  are  designed  to 

'/.    run  using  variables  referenced  to  an  r-sub-D  (ro)  which  is  equal  to  the 

'/.    size-i  Li  of  the  phase  screeni  i-e-i  D  is  actually  incorporated  as  D/ro- 

V. 

'/.    This  script  will  utilize  the  same  calculations  and  fix  the  aperture 

'/.    sizei  Di  to  a  fixed  sub-aperture  sizei  di  and  run  through  a  range  of 

'/.    (ro)  values-  The  'fixed'  parameters!  Di  di  and  Li  therefore  must 

'/.   be  scaled  for  each  run  through  different  (ro)  values- 

y. 

'/.    The  scintillation  has  been  incorporated  alsoi  to  run  through  a  range  of 

'/.    values  using  sigma-sub-L  squared- 

'/. 

y. 

'/.  1)  N  =  lb  thru  IDEM  possible  choices  (in  powers  of  2) 

'/. 

'/.  2)  Required  Matlab  functions  or  scripts  for  this  to  work  at  compile  time 

'/.  -  -  Scinfun-m 

'/.  --    avint-m 

'/.  —  hv57-m 

'/.  --psdgen-m 

'/. 

'/.  3 )  Choose  desired  values  for: 

'/.  --    ro  Erange  of  valuesHi 

'/.  --  sigma  Erange  of  valuesli 

'/.  --  Di  aperture  diameteri  in  metersi 

'/.  --  di  sub-apertures  effective  diameteri  in  meters1. 

'/.  --  Li  effective  length  of  phase  screeni  in  metersi 

'/.  --  fni  number  of  loops  for  desired  computation 

'/.  --  filenamei  for  saving  data 

V. 

y.y.y.  */.  y.  y.  */.  y.  y.  v.'/.'/.'/.  y.  y.  ■/.  y.y.y.y.  v.  y.  v.  '/.v.'/.'/.v.'/.'/.'/.'/.'/.'/.'/.  y. '/.  y.  y.  y.  y.  y. '/.  y.  y. '/.'/.'/.'/.  y.  y.  y.  y.  z '/.  y.  y.  y.  y.  '/.  y.  y.  y.  y.  '<  y.  y.  y.  v.  >.  y.  y. 

tD  =  clocki 

'/.    Define  parameters  for  Ro  manipulation 

N  =  E5bi 

ro  =  C-DM0i.D37i.D3mi 
sigma  =  E ■ 11 i -D7i • 03i -Dili 

D  =  1-li  '/.  aperture  diameteri  in  meters 

d  =  D.li  y.    sub-aperture  diameteri  in  meters 

L  =  3*i  y.    screen  lengthi  in  meters 

aa  =  length (ro) i 

bb  =  length ( sigma) \ 

'K    Scale  the  physical  parameters  to  the  varying  Ro  value 
LL  =  L./roi 
DD  =  D./roi 
dd  =  d • /roi 
clear  L  D  di 
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'/.    Calculated    results    for    S-functions    in    denominator    of    alpha 

Sl  =  EE.E7mflfll4be-07-.    3  •  MTlS^^Se-DT-,    5  .  M333D535e-ll-,    8  .  14808817^6-13    ... 

1.35H7,!Hflc1e-m-.    E-DbTfibTfibe-lb-,    3.S34117DDe-ia3i 
SE  =  EE.EL.3M87^Ee-05-,    5  .  SSTTTlMbe-D?-.    1 .  153bD5M3e-08-,    E  •  MblSfi3ne-10    ... 

S.m3S1flHbe-12i    l-D5TDfi7Dbe-13-.    E .  157bD7clfie-lSl!--, 
sl=Sl(logE(N)-3)^ 
sE=SE(logE(N)-3)n 
clear    SI    SE't 

M    =    N/Ei 

alpha  =  sqrt  (b.  853677*  ((E*I1)  a  (-5/3) -D- =18=1^55*  (E*M)  a  (-E)  )  /  (  sE-sl)  )  =, 

'/.    Define  the  spatial  frequency  range  for  phase  screen  generation  later 

f  =  CD:  (M-l)  C-I1):(-I)3i 

fE=f -AEi 

FE  =  ones(E*M-.l)*fEi 

FE=FE+FE' \ 

FE(l-.l)=eps:> 

clear  f  fSi 

F3  =  alpha  *  (FE • a (-11/1E)  )  i 

F3(l-il)=0n 

clear  alpha 

PpUt=Cli 

IntenRef =C3i 

for  qq  =  1 :  aa 

'/.    Generate  the  spatial  frequency  scaled  by  1/LL  and  then  in  E-D 

kappa  =  CD:((N/S)-1)  (-N/E) : (-1) 3/LL ( qq) i 

kappa  =  meshgrid (kappa) t 

kappamag  =  sqrt (kappa- AE+kappa '• aE) \ 

kappa  =  1 inspace (min (min (kappamag )) Tmax (max (kappamag )) nE*N ) i 

'/.    Generate  the  spatial  freq.  transform  of  the  aperturesi  which  again  has 

'/.    been  scaled  by  the  varying  Ro  value 

x  =  pi*kappa*DD (qq) \ 

lilt  =  D.S5*pi*(DD(qq)AE)*(Cl  E*bessel  j  ( 1  -.x  (E  :  (E*N)  )  )  .  /x  (E  :  (  E*N)  )  1)  i 

'/.    Calculate  the  reference  intensities  for  Strehl  computation  later 

ppliltQ  =  spline  (kappa-ildt )  i 

U  =  sqrt  (E)  *ones  (N)  \  '/.    Phase  screeni  U 

Ut  =  fftE(U)^  '/.    Phase  screen  transform-!  Ut 

lilt  =  ppval  (ppliltOikappamag)  \       '/.    Aperture  transform!  Ut 

V  =  abs  (  if  f  tE  (Ut  •  *Ut)  )  i       '/.    Uniform  reference  intensity  for  diff-limit 

"/^Concatenate  results  for  later  multiple  use 
IntenRef  =  [IntenRef-.  V  (1-.1)  aE]'-, 
ppblt  =  Epplilti  ppldtO  3i 

end 

clear  x  lilt  ppUtD  U  Ut  V  qqi 

Prob  =  sparse(l->l-.D-iaa*bb-.10e3-il0e3*aa*bb)i 
uu  =  1 :  leM  =. 

randn( 'state ' iSum(10D*clock) ) i 
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'/.   Begin  major  looping  through  ranges  of  (ro)  and  (sigma)  values 
for  fn  =  l:750i 


count  =  li 
for  nn  =  1 : aai 
for  mm  =  l:bbi 


'/.    Keep  track  through  the  interior  loops 
'/.    Loop  through  the  (ro)  values 
'/.    Loop  through  the  (sigma)  values 


scrn-fn  =  randn (E*M  iE*M ) +i*randn (E*M  iE*f1 ) i 
scrn— fn  =  F3 • *scrn— f ni 

x=pi*sqrt(FE)*dd(nn)/LL(nn); 

x  ( l-il )  =eps  i 

T  =  sqrt(l-(£*bessel  j(l-.x)  ./x)  .  aE-  (  M*bessel  j  (2,x  )  • /x)  .  »2)i 

clear  xi 

T(lil)=0i 

scrn  =  scrn— fn.*Ti  clear  Ti 

scrn  =  ifftE(scrn)  *  ( (LL (nn ) ) a (5/b ) ) \ 

'/.    Function  call  to  get  the  scintillation  screen  scaled  to  sigma  value 
CscrnScl  =  Scinf un (N nsigma (mm) ) \ 

kappa  =  CD:((N/E)-1)  (-N/H) : (-l)3/LL(nn) \ 

kappa  =  meshgrid (kappa) \ 

kappamag  =  sqrt (kappa • AE+kappa '• aE) \ 

U  =  sqrt  (E)  *(exp(  i*real  (scrn) ).  *exp(  (real  ( scrnSc) -sigma  (mm)  )))  =. 

Ut  =  fft2(U)i 

lilt    =    ppval  (pplilt  (  nn  t  :) -i  kappamag  )  i 

V  =    abs(ifftE(Ut.*Wt) ) i 

prob  =  round(10e3*(V( : ) . AE/(IntenRef (nn) ) ) ) i 

ii  =  find(prob  ==  0)"-. 

prob ( i  i )  =  li 

Prob  =  Prob  +  sparse  (count  iprob  -i]naa*bb-il0e3)  \ 

U  =  sqrt(E)*(exp(i*imag(scrn) ) . *exp ( (real(scrnSc)-sigma(mm) ) ) ) i 
Ut  =  fft£(U)i 

V  =  abs(ifftE(Ut.*li)t)  )  \ 

prob  =  round(lDe3*(V( : ) . AE/(IntenRef (nn) ) ) ) i 

ii  =  find(prob  ==  0)=. 

prob(ii)  =  1 \ 

Prob  =  Prob  +  sparse  (count  iprobil  -iaa*bb-il0e3)  \ 

count  =  count+li 

end 
end 

'/.   Save  the  data  every  few  runs  in  case  of  unsolicited  stopping- 
if  mod(fn-iED)=  =  D 

ttt=etime(clock-itO)i 

save  Sin7Dl  Prob  fn  ttti 

disp(  C  'frame  *'  intSstr(fn)  '.'  1    )'■, 
end 


end 

clear  nn  scrn-fn  scrn  FE  F3  ppldt  kappamag: 

clear  U  Ut  lilt  V  kappa  prob  scrnSci 
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fn  =  E'done  '  intEstr(fn)  '  runs'J 

ttt=etime(clockitO) i 
uu  =  uu/le4i 

'/.    Save  the  data  here  at  the  end 
save  Sin701  Prob  fn  ttt  uui 

'/.    Execute  TimeDisplay  M-file-i  if  it  exists- 
if  exist (' TimeDisplay • m ' ) 

TimeDisplay(ttt) i 
end 

'/.    Important  to  have  this  'quit'  when  running  on 

'/.    Silicon  Graphics  SGI  machines 

quit 
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APPENDIX  B.    REQUIRED  MATLAB  FUNCTIONS  FOR  SIN.M' 


Scinfun.m 

'/. 

'/.    Function!  called  Scinfun-.  to  incorporate  Dr.  Fried's  two  programs-. 

'/.    psdgen  and  log— amp-screen  •  As  the  adaptive  optics  are    phase-correcting 

'/.  only-,  no  corrections  will  be  introduced  here-   The  turbulence  will  be 

'/.  varied  through  the  severity  of  the  C-sub-N  squared!  which  in  turn  will 

'/.    be  set  to  the  values  used  in  the  Turro  script  runs-  This  will  also 

'/.    require  a  wavelength  selection!  not  needed  in  Turro- 

'/. 

y. 

'/.    N  =  It  thru  IDEM  possible  choices  (in  powers  of  E) 
y. 

function  Cscrn  n  siglsql  =  Scinf  un  (N  ->ro  ) 

'/.    Define  parameters  for  R  o  manipulation 

lamb  =  0-fiMe-bi      'A    AMD  nanometers  -  near  infra-red  light  propagation 

D  =  D-Mi  y.    physical  aperture  diameter  (m) 

L  =  3i  '/.    physical  phase  screen  diameter  (m) 

'/.    Scale  the  physical  parameters  to  the  varying  Ro  value 

LL  =  L/roi  '/.    all  arrays  (1  x  length(ro)) 

DD  =  D/roi 

clear  D  i  '/.    keep  (ro)  for  later  conversion  to  actual  value  (cm) 

'/.  Generate  the  PSD  for  the  scintillation  in  the  turbulence 
air  =  logspace(D-.loglO(3eM)  -.ESQ)  i 
CNEO  =  hv57(air)i 

'/.  Still  need  to  coordinate  the  ro  with  CNE  --  basically  dial  the 
'A    CNE  value  using  the  ro  value  as  the  knob 

'/.    bie  will  assume  that  (ro)  from  hvS?  equals  5-7 
hveeS?  =  -05?:. 

randnC  'state'  -i sum (10D*clock)  )i 

aaa  =  ( (ro*L) /hvee57) * (-5/3) ; 
CNE  =  aaa    *    CNEDn 

Csiglsq-.ppPSDH  =  psdgen  (air-.  CNE  t  lamb  )n 

'/.    Define  the  spatial  frequency  range  for  phase  screen 

mt  =  1:  (N/E)  i 

mt  =  Cmt-1  (  (  (N/E)+l)  :N)-N-13i 

kappa  =  mt/LLi 

kappa  =  meshgr id (kappa) i 

kappamag  =  sqrt  (  kappa  •  *£  +  kappa'  .aE)=. 

clear  kappa '-. 
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srPSD  =  sqrt (ppval (ppPSD t kappamag ) ) i 
clear  kappamagi 

alfa  =  NaE/LL. 

scrn— fn  =  randn ( N ) +i*randn  (  N  )  i 
scrn— fn  =  (alf a*srPSD) ■ *scrn— f ni 
clear  srPSDi 

scrn  =  if  f  tE  (  scrn— f  n  )  =. 


avint  .m 

'/. 

'/.    Function!  called  avinttxiynxloixup)  to  perform  integration 

'/.    of  a  curve  defined  by  data  points  (xiy)  using  overlapping 

'/.    parabolas-  Integrates  from  xlo  to  xup-  These  do  not  need  to  be  one  of 

'/.    the  actual  abscissa  pointsi  x(i)-  Must  have  xlo  <  xup-  These  default 

'/.    to  xlo  =  x(l)  and  xup  =  x  ( length  (x ))  i 

'/. 

'/.    Coded  from  Davis  &  Rabinowitzi  "Methods  of  Numerical  Integration!" 

'/.    pg.HA3  --  *tjb  fi/1/11 

'/. 

'/. '/. '/. '/. '/. '/. '/. '/. '/. '/.  y.'/.'/.y.'/.  '/.  y. '/.  v.'/.'/.v. '/.  y.  y.  y.  y.  y.  y.  y.  y.  v.  y.  y.  y.  y.  y. '/.  y.  y.  y.  v.  y.  y.  y.  v.  y.  v.  y.  y.  y.  y.  y.  y.  y.  y.  y.  y.  y.  y.  y.  y.  y.  y. '/.  y.  y.  y.  y. '/.  y. '/. 

function  I  =  avint  (  x  -,y  -,  xloixup  ) 

n=length ( x  )  ^ 

if  narginOn  xlo  =  x  (1 )  ixup  =  x  (  n  )  i  end 

syl  =  xlo  =. 

for  ib  =  l :  n 

if  x(ib)<=xloi  breaki  end 
end 
ib  =  min(max(E->ib)  -in-1)  i 

for  j  =  n : -1 : 1 

if  xup>=x(j)i  breaki  end 
end 
j  =  max(ib-imin(j-in-l)-l)T 

sum  =  D:i 

for  jm=ib :  j 

xl=x(jm-l) ix2=x(jm) ix3=x(jm+l) i 

tl=y( jm-l)/( Cxl-x5)*(xl-x3)); 

tH=y( jm)/( (xE-xl)*(xE-x3) ) i 

t3=y( jm+l)/( (x3-xl)*(x3-xE) ) i 

A=tl+tE+t3i 

B=-(xE+x3)*tl  -  Cxl+x3)*tE  -  (xl+xE)*t3i 

C=xE*x3*tl  +  xl*x3*tE  +  xl*xE*t3; 

if  jm==ib 

ca  =  Aicb  =  B:iCC  =  Ci 
else 

ca=.5*(A+ca)icb=-5*(B+cb)icc=.5*<C+-cc)  i 
end 
sum  =  sum  +  ca*  (  xEA3-sy  1  *3  )  /3  +  cb* • 5* (xEAE-sy IaE )  +  cc*  (  xE-syl ) "-. 
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ca  =  A;  cb  =  B  =i  cc  =  Ci 
s  y  1  =  x  E  i 


end 


I  =  sum  +  ca*(xupA3-sylA3)/3  +  cb* . 5* (xupAE-syl aE)  +  cc* ( xup-sy 1 ) i 


hv57  .m 
v. 

'/.    Function!  called  hv57i  computes  the  HV5/7  CnAE  profile  -  This  was 

'/.    converted  to  m  from  the  Fortran  source  in  the  tOSC  library. 

'/. 

'/.    alt  is  an  array  of  values  and  the  function  returns  an  array  of  values- 
'/. 

'/.  MODEL  =  XXXD7 

'/. 

'/.  modelname  =  'Modified  Hufnagel-Valley' 

V. 

function  cnE  =  hv57(alt) 

model  =  E107=i 

cnE  =  0*alti 

inatmos  =  (0  <=  alt  <=  £4000):. 

if  all (-inatmos ) 

return 
end 

al  =  al t ( inatmos ) i 

v  =  ceil (model/100) \ 

cnE(inatmos)  =  (S • Ee-S3*al . a10  •*  ( ( v/E7 ) AE*exp ( -al/1000 ) )  ... 

+    l.e-lb*exp(-al/1500) )*2-7i 
low    =    (0    <=    alt    <    11000) i 

cnE(low)     =    cnE(low)+1.7e-14*exp(-alt  (low)/lD0)  =. 


psdgen . m 

'/. 

'/.  Function!  called  psdgen->  used  to  generate  results  for  the  PSD  of  the 

'/.  log-amplitude  covariancei  Phi-sub-1  ( k )  i  and  for  the  log-amplitude 

'/.  variance!  s  igma-sub-1-squared  •  It  takes  as  inputs  values  for  the  wave 

'/.  lengths  and  the  value  of  C-sub-N-squared  at  various  positions!  S! 

'A  along  the  propagation  path-  The  PSD  result  is  reported  out  in  spline 

'/.  formi  suitable  for  use  by  the  function  ppval  • 

'/. 

y.  y.  y.  y.  y.  v.  y.  y.  y.  y.  y. '/.  y.  y.y.y.  y.  y.  y.  y.  y.  y.  y.  y.  y.y.y.y.  y.  y.  y.  y.  y. '/.  >.  y.  z  y.  y.  y.  y.  y.  y.  v.  y.  y.  y.  y.  >.  y.  y.  %  y.  y.  '/.'/.'/.v.  y.  y.  y.  '<  y.  y.  >.  v.  y.  '<  y.  y.  y.  y. 


63 


INPUTS 
s        dim) 
CNE      dim) 

lambda   dil) 

OUTPUTS 
siglsq   dil) 

ppPSD    diSm) 


Position  along  propagation  path  (in  meters)- 
Value  of  C-sub-N-squared  at  positons  along  the 
propagation  path  (in  mA(-E/3))- 
Uave  length  (in  meters) • 


Log-amplitude  variancei  sigma-sub-1-squared 

( in  NpAE ) • 

Spline  coefficients  for  the  evaluation  of  the 

power  spectral  density at  a  spatial  frequency 

whose  magnitude  is  given  in  cycles/meter 

(in  NpAE/ ( cycles/m) a2 ) . 


7. 
7. 
7. 
7. 
7. 
7. 
7. 
7. 
7. 
7, 
7. 
7, 
7. 
7. 
7. 

function  EsiglsqippPSDJ=psdgen(siCNEi lambda) 

k  =  2*pi/lambdai 
F  =  'avint' i 

B  =  EA(-t4/3)*piA(3/2)/(3*cos(pi/lE)*gamma(E/3))i 
siglsq  =  B*k*(7/b)*  f eval (Fis iCNE . *s  .  a (5/b ) ) i 

A  =  (EA(5/3)*sqrt(pi)/T)*(gamma(ll/b)/gammadD/b)  )  i 

kappamax  =  SO*  (  k/ (E*pi*max  (s  )  )  )  =. 

kappa  =  kappamax*logspace  (-5iDi333) i 

kappamll3  =  A*kAE* ( E*pi ) a (-fl/3 ) *kappa • a (-11/3) i 

kappasq  =  (E*pi ) AE*kappa . AE/k \ 

PSD  =  zeros(li333)i 
for  n=l:333 

PSD(n)=  kappamll3(n)*feval(FiSiCNE.*d-cos(kappasq(n)*s)))i 
end 

ppPSD  =  spl ine (kappaiPSD) i 
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APPENDIX  C.  EXTRA  INTERESTING  RESULTS 

Along  the  course  of  completing  this  thesis,  not  everything  went  according  to  plan, 
nor  as  one  would  want  to  happen.  The  following  mistake  was  the  cause  of  nearly  a 
month  in  lost  time,  but  hopefully  at  a  cost  of  deeper  understanding.  At  some  point  during 
the  preparation  phase  of  the  combined  phase  and  scintillation  program,  "sin.m,"  the 
author  inadvertently  input  unrealistic  values  for  the  aperture  diameter  and  its  sub-aperture 
diameter.  The  "unnoticed  until  examination  of  results"  mistake  led  to  an  unrealistic  D/d 
ratio  of  4.0  and  unintentionally  broadened  the  Strehl  range  of  the  Gaussian-shaped  bell 
curves.  The  resultant  "straight"  cumulative  probability  lines  were  noticeably  not  straight, 
and,  in  fact,  the  "sin.m"  program  could  not  duplicate  the  previous  "turro.m"  results  when 
zero  scintillation  was  initialized.  After  a  line-by-line  check  of  the  two  programs,  the 
different  aperture  initializations  were  noticed,  but  not  before  one  billion  simulation 
results  were  developed  and  examined.  The  results  from  that  data  set,  called  at  the  time 
"combine.m,"  are  reproduced  here  in  the  appendix  for  the  perusal  of  interested  readers. 

As  presented  in  the  body  of  the  thesis,  the  following  Table  C-l  is  a  concise 
summary  of  the  initialized  values  for  this  set  of  results.  It  can  be  seen  that  we  ran  five 
values  of  r0  versus  seven  total  values  of  g(  (including  the  value  of  zero  scintillation). 

The  following  Figures  (C1-C3)  should  be  self-explanatory  after  the  thesis 
presentation  for  the  subsequent  runs  of  'Sin,"  although  we  will  mention  two  points  of 
interest.    First,  the  range  of  Strehl  ratio  realizations  is  quite  large,  from  approximately 
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0.25  up  to  0.75.    This  can  be  physically  explained  as  having  effectively  only  four  sub- 
apertures 


Parameter 

Initialized  Value 

Aperture  Diameter 

0.40  (m) 

Sub-Aperture  Diameter 

0.10  (m) 

D/d  Ratio 

4.0 

Screen  Diameter 

3.0  (m) 

Wavelength  (X) 

0.840  (microns) 

Values  ofr0  (in  meters) 

0.1260,  0.1200,  0.1 140,  0.1080,  0.1020 

Resultant  aVr0  Ratios 

0.7937,  0.8333,  0.8772,  0.9259,  0.9804 

Values  of  Scintillation  (in  Np  ) 

0.10,  0.08,  0.06,  0.03,  0.02,  0.015,  0.00 

Screen  Size(N) 

256  x  256 

Table  C-l.  Initialized  Parameters  for  Results  of  Combine'  Computer  Runs,  10-26  September  1998. 

across  the  entire  aperture,  not  able  to  provide  much  correction  at  all.  Second,  as  a 
consequence  of  such  a  wide  spread  in  Strehl  ratio  realizations,  the  probability  lines  in 
Figure  (C2)  are  very  long  and  not  as  straight  as  noticed  in  all  of  the  other  runs  (which 
were  done  with  a  more  realistic  D/d  ratio  in  every  case).  This  can  be  physically  explained 
as  approaching  the  linear  boundaries  of  zero  and  one  on  either  side  of  the  Strehl  axis. 
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Figure  CI.  10-26  September  1998  'Combine'  Computer  Code  Results 
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Figure  C2.  Strehl  vs.  Probability  Lines  at  d/r0  =  0.794 
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Figure  C3.  Strehl  vs.  Probability  Lines  for  Ge  =  0.06  Np2 
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